Non-uniform liquid— crystalline phases of parallel hard rod-shaped particles: From 

ellipsoids to cylinders. 

Y. Martinez-RatorQ 

Grupo Interdisciplinar de Sistemas Complejos ( GISC), Departamento de Matemdticas, 
Escuela Politecnica Superior, Universidad Carlos III de Madrid, 
Avenida de la Universidad 30, E-28911, Leganes, Madrid, Spam 

E. VelasccQ 

Departamento de Ftsica Teorica de la Materia Condensada and Instituto de Ciencia de Materiales Nicolas Cabrera, 

Universidad Autonoma de Madrid, E-28049 Madrid, Spain 
(Dated: February 26, 2008) 

In this article we consider systems of parallel hard superellipsoids, which can be viewed as a possi- 
ble interpolation between ellipsoids of revolution and cylinders. Superellipsoids are characterized by 
an aspect ratio and an exponent a (shape parameter) which takes care of the geometry, with a = 1 
corresponding to ellipsoids of revolution, while a = oo is the limit of cylinders. It is well known that, 
while hard parallel cylinders exhibit nematic, smectic, and solid phases, hard parallel ellipsoids do 
not stabilize the smectic phase, the nematic phase transforming directly into a solid as density is 
increased. We use computer simulation to find evidence that for a > a c , where a c is a critical value 
which the simulations estimate to be in the interval 1.2-1.3, the smectic phase is stabilized. This is 
surprisingly close to the ellipsoidal case. In addition, we use a density-functional approach, based on 
the Parsons-Lee approximation, to describe smectic and columnar ordering. In combination with a 
free-volume theory for the crystalline phase, a theoretical phase diagram is predicted. While some 
qualitative features, such as the enhancement of smectic stability for increasing a, and the probable 
absence of a stable columnar phase, are correct, the precise location of coexistence densities are 
quantitatively incorrect. 
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I. INTRODUCTION 

Hard interaction models have played an important role 
in the understanding of the nature and structure of sim- 
ple liquids and crystals made of particles with spherical 
symmetry. For anisotropic particles, the hard ellipsoid 
(HE) and hard-spherocylinder (HSC) models have played 
a role similar to that of the hard-sphere model (HS) in 
simple liquids, though these models are not so univer- 
sal as the HS model. In particular, HE and HSC fluids 
exhibit an isotropic-nematic phase transition but, while 
the HSC fluid shows a stable smectic phase, all evidence 
to date suggests that the HE fluid does noli Clearly, 
the formation of the smectic phase must be the result of 
delicate packing effects directly related to particle shape. 
This feature of the HE fluid is a clear disadvantage in the 
formulation of perturbation theories for liquid crystals^. 

The problem of why hard ellipsoids do not get sta- 
bilised into a layered smectic structure is in intriguing 
one. The properties of a perfectly aligned fluid of ellip- 
soids can be mapped onto those of hard spheres, and the 
HS fluid does not exhibit phases with order intermediate 
between the fluid and the crystal. On top of that, orienta- 
tional freedom probably plays against smectic formation. 
Recent simulation work on parallel hard ellipsoids aug- 
mented by an isotropic square well have shown that, in 
this system, smectic layers can be stabilised-; however, in 
view of the rather artificial model potential used, the re- 
sult probably does not reflect any essentially interesting 



underlying property of ellipsoids. 

The physical reason for the absence of smectic order in 
the HE fluid obviously lies in the geometrical properties 
of an ellipsoid. Wen and Meyer— addressed the general 
problem of smectic formation in fluids of parallel hard 
rods. They provided an explanation in terms of the en- 
tropy gain involved in the increased packing efficiency in 
smectic layers, which more than outweighs the entropy 
loss associated with the onset of layering with respect 
to the nematic phase. This efficiency is very much re- 
duced in ellipsoids due to particle shape, since ellipsoids 
arranged in a layer leave too much void space. Alter- 
natively one may think that increasing packing, which 
would involve filling this space, entails interlocking be- 
tween the layers, which promotes crystalline order but 
discourages formation of the smectic phase. 

A question that can be asked to understand this prob- 
lem from a different perspective is the following: if we 
perturb the shape of an ellipsoid in the direction of a 
spherocylinder or a cylinder, both of which exhibit smec- 
tic phases^, for which particle shape does smectic sta- 
bility set in? The answer to this question may provide 
some further insight into the relation between particle 
geometry and smectic stability. As a bonus, it would 
help formulate more useful hard-body models that can 
be used in perturbation theories. 

Only a few previous studies have addressed this is- 
sue. Of particular relevance to our study is that of 
Evans&, who used an Onsager second-virial coefficient 
approximation to investigate smectic formation in flu- 



2 



ids made of hard ellipsoids, hard spherocylinders and 
hard ellipocylinders, a particle with a shape somehow 
intermediate between that of the first two. The ellip- 
soids, both parallel and with unconstrained orientations, 
did not form a smectic phase before the crystal, while 
the other particles did at some particular density prior 
to crystallization. It was concluded that ellipsoids are 
pathological in that they do not form a smectic phased. 



In the present paper we again address this problem by 
studying a continuum of particle shapes, in the limit of 
parallel particles, but this time interpolating between the 
ellipsoid and the cylinder by means of a model, the hard 
superellipsoid (HSE) of revolution, containing an expo- 
nent a that can be varied continuously. Monte Carlo 
(MC) simulations are used to analyse the stability of the 
smectic phase and other phases with partial positional or- 
der (i.e. columnar phase) of the parallel model (PHSE), 
in the region of geometries close to the ellipsoid, and 
an approximate phase diagram as a function of a is ob- 
tained. The columnar phase is not stable in the phase 
diagram. MC simulations of freely rotating cut spheres 
have shown that the columnar phase can be stabilised 
for some range of aspect ratios^. However simulations of 
parallel cylinders do not give conclusive evidence for the 
presence of columnar symmetry in the phase diagram 5 , 
and a recent fundamental measure theory (FMT) for par- 
allel cylinders 8 also rules out the columnar phase as a 
thermodynamically stable phase. However, one of the 
conclusions of the present paper is that the smectic phase 
can be stabilised with respect to the nematic phase for 
some type of PHSE particles. Therefore, in the approxi- 
mation of parallel particles, ellipsoids do not seem to be 
pathological; rather it is a continuum of particle shapes, 
close to the ellipsoidal, that do not exhibit smectic order. 
To complement these studies, a density-functional the- 
ory based on the Parsons-Lee approach^ was used to in- 
vestigate the relative stability of the smectic phase with 
respect to other phases with lower symmetry. For the 
crystal phase, a free-volume approximation is used. The 
main conclusion that can be drawn is that the Parsons- 
Lee approximation, in combination with free-volume the- 
ory, only gives a qualitative picture of the relative stabil- 
ity of the phases. The precise location of the transition 
points are incorrect. Other features, such as the fact that 
the columnar phase is not stable, appear to be correctly 
described. Resort to other, more sophisticated theories, 
e.g. of the FMT type, would be needed to obtain a bet- 
ter description. However, although the initial framework 
of the FMT for general convex particles was initially de- 
veloped by Rosenfeld^, its application only works for 
isotropic fluids, as it was shown in Ref.— . For some sim- 
ple geometries (such as parallelepipedic) the FMT can 
be formulated from first principles with restricted par- 
ticle orientations^. The generalisation of the FMT to 
geometries like the ones described here seems to be a 
rather difficult task. 



II. PARTICLE MODEL 

We consider parallel particles with symmetry of revo- 
lution (axially symmetric about the z axis). The superel- 
lipsoid results by making a superellipse in the x-z plane 
revolve around the z axis, with the equation 
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a and b are the semi-lengths of the particle in the xy 
plane and along the z axis, respectively, and a is the 
shape exponent. This equation defines a particle with 
a geometry interpolating between the ellipsoid (a = 1) 
and the cylinder (a = oo). Note that the PHSE model is 
different from that used by Evans&, whose model inter- 
polates between an ellipsoid and a spherocylinder. Fig. 
Q] depicts a few bodies (with a = b, which could be called 
superspheres) for different values of a. Thermodynamic 
and structural properties of such a fluid of parallel par- 
ticles scale trivially with particle elongation, so that it is 
sufficient to consider the case a = b = <r /2, with cr the 
particle breadth. 







Figure 1: Superellipsoids with a — b (superspheres) for values 
of a equal to 1, 1.5, 2 and 5 (from left to right). 

The stable low-density phase of the PHSE model is a 
nematic phase (all particles oriented along the z axis and 
with their centers of mass disordered positionally). At 
high densities we expect a solid phase with some crys- 
talline structure (the question of which structure is sta- 
bilised will be addressed later). At intermediate densities 
we could, in principle, expect different phases with partial 
order (smectic or columnar) to get stabilised, depending 
on the value of a. 



III. THEORETICAL TOOLS 

To explore the phase behaviour of the PHSE model we 
used constant-pressure MC simulation. We have simu- 
lated systems of N = 1-2 x 10 3 particles, with a = 1.0 
2.0. Both compression runs from the low-density nematic 
phase and expansion runs from a high-density crystalline 
structure were conducted. From these, equations of state 
are obtained, and the different phases may be identified. 
Relative stability between nematic, smectic and crystal 
phases (some of which, as we will see, undergo a first- 
order phase transition) cannot be ascertained, since no 
attempt was made at computing absolute free energies. 
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Structural quantities, such as density profiles along the 
ordering direction, p(z), and radial distribution functions 
parallel and perpendicular to this direction, g\\(z) and 
g±(R), were also calculated. All this information allows 
us to get a picture of the trends in phase behaviour that 
can be expected as the exponent a is varied. 

To complement simulation results, a number of the- 
oretical analyses have been carried out. In the high- 
density region, the classical free-volume theory has been 
implemented. The free- volume free-energy density is 
written as 
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where N is the number of particles, i> free the free volume 
available to a particle when the others are kept fixed in 
their lattice positions, (3 — 1/kT and A is the thermal 
wavelength. For parallel hard cylinders this is an an- 
alytical function of density. In the general case this is 
probably not true, and we have computed w free by MC 
integration. In the low-density region, a simple virial 
expansion for the excess free energy, 
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has been considered to calculate the nematic-smectic 
bifurcation lineA 3 -, with an aim to comparing with the 
Parson-Lee approach (see later) . In the above expansion 
/(r) is the Mayer function of two parallel PHSE, and p{v) 
the local density distribution. This expansion is meaning- 
ful only for the low-density nematic phase and expected 
to rapidly fail as the system density is increased, but at 
least it may give an indication as to whether the system 
is prone to developing smectic ordering and, if so, how 
this tendency depends on particle shape. Specifically, we 
have applied a bifurcation analysis based on the above 
expansion (presented in Appendix A), using second- and 
third-order terms in p(r). At higher density, i.e. in re- 
gions where the smectic or columnar phases may be sta- 
ble, an alternative is to use a Parsons-Lee (PL) scheme^. 
In the PL approach, an approximate resummation of the 
virial series is performed, using the exact second virial 
coefficient. That is, we take 
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where B^ s are hard-sphere virial coefficients. With this 
scaling of the PHSE virial coefficients, the excess free 
energy writes 
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where 'J'hsW is the excess free energy per particle of 
a HS fluid of the same packing fraction r\ as our fluid 



of PHSE particles. The packing fraction is given by 
V = PqVq, with po the mean density and vq the parti- 
cle volume, given by 
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where B(x, y) is the beta function. For vE , hs( ? 7) we use the 
Carnahan-Starling expression, 5'hs( ? 7) = (4 — 377)77/ ( 1 — 
77) . In the minimisations of the total free energy, 

m = r 1 J drp(r) {log [p(r)A 3 ] - 1} 

+ ^ex[p], (7) 

where •7 r e x[/>] is given by ((3]) or (|5]) , the density profile p(r) 
is parametrised in some convenient way. For the smectic 
and columnar phases at low densities we use a truncated 
Fourier expansion. For high densities a Gaussian param- 
eterization is used. Details of these calculations are given 
in Appendix B. 

IV. RESULTS 

In this section we use the packing fraction 77 as a con- 
venient measure of density. For a given density, since v 
is an increasing function of a, r\ increases slightly from 
the ellipsoid to the cylinder. In Fig. [21 MC data for 
reduced presure pvo/kT versus packing fraction 77, for 
all the particle shapes considered, are shown (from now 
on we only show simulation results for the systems with 
N ~ 10 3 particles; selected checks with twice as many 
particles did not give any quantitative differences) . In all 
cases the low-density phase is a nematic (N phase), since 
all particles are parallel, whereas the system crystallises 
at high density (K phase). We begin by discussing the 
high-density region, where the crystal phase is the stable 
phase. 

As expected, compression runs starting from the low- 
density nematic phase (open squares in Fig. [2j) ultimately 
create a defected crystal, via a first-order phase tran- 
sition (this appears as a clear density discontinuity in 
all cases). Expansion runs from a perfect high-density 
crystal (filled squares in Fig. [2]) give a solid branch 
with defect-free structures in most instances. However, 
in small systems defects consisting of columns of parti- 
cles that have moved along the director by a fraction of 
the unit cell are created. This may give the impression 
that averaged structures correspond to a stable columnar 
phase. A related situation was found in the simulations of 
Veerman and Frenkel on parallel hard cylinders^, where 
a strong dependence of columnar stability on system size 
was found. 

To investigate this, we prepared initial configurations 
with columnar symmetry at high density; in expansion 
runs, the system always stayed as a columnar phase (tri- 
angles in Fig. [21. Therefore, for large system sizes, 
systems with columnar and crystalline structures look 
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Figure 2: Equation of state pressure-packing fraction for var- 
ious values of shape parameter a (shown in each panel), as 
obtained from constant-pressure Monte Carlo simulations. 
Open squares: compression runs starting from nematic phase. 
Filled squares: expansion runs starting from crystalline phase 
with ABC symmetry. Triangles: expansion runs starting from 
columnar symmetry. Identification of phases is made with 
symbols (N, S and K for nematic, smectic and crystal, re- 
spectively). 
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Figure 3: Radial distribution functions along (continuous 
line) and perpendicular to (dashed line) the director, g\\ (z) 
and g±(R), for the PHSE fluid with a = 1.5; (a) 7} = 0.38; 
(b) T) = 0.43; and (c) 77 = 0.53. 



as though they are actually mechanically stable, and 
there seems to be a large free-energy barrier between 
the two structures. However, compressions from the ne- 
matic never give rise to configurations containing par- 
tial columnar order (which would seem easier to gener- 
ate than full three-dimensional order). Without explicit 
calculation of free energies, no definite conclusion can be 
reached on the relative stabilities of columnar and crys- 
talline phases; however, given that no clear evidence for 
columnar ordering has been found in simulations of par- 
allel hard cylinders^, and that FMT calculations on the 
same system predict that columnar order is not thermo- 
dynamically stable^, we believe it unlikely that a colum- 
nar phase may get stabilised in PHSE. 

As the crystal is expanded to lower density the sys- 
tem looses translational order (either totally or partially), 
and becomes fluid via a discontinuity in density. Fig. 
[3^c) shows the radial distribution functions along and 
perpendicular to the director, gn(z) and g±(R), in the 
crystal phase for the case a = 1.5. Both present a high 
degree of structure, as expected in a crystalline phase. 
Therefore all evidence suggests that there is a first-order 
fluid-crystal (i.e. freezing) transition. 

The question on the nature of the crystalline phase 
and how the system is prepared in the expansion runs 
deserves some comments. In fact the symmetry of the 
structure seems to change at some value of a. The generic 
crystal structure consists of stacked triangular layers, ei- 
ther in phase (AAA structure) or out of phase. The latter 



may have ABCA..., ABAB..., or random-stacking struc- 
tures, all of these being compact structures (i.e. they 
share the same value of the close-packing density) . A pri- 
ori one could consider the AAA, ABC and ABAB struc- 
tures to be the most stable candidates, reflecting simple 
hexagonal-, face-centred-cubic- and hexagonal-close- 
packed-like symmetries, respectively (the final structures 
obtained by the simulations do not have these exact sym- 
metries, since the intralayer unit-cell distance does not 
exactly correspond to the interlayer spacing expected in 
these structures; obviously this is a consequence of the 
asymmetry of the particles along and perpendicular to 
their symmetry axes). The situation, for lack of a more 
detailed analysis based on free-energy estimations is, as 
usual, uncertain (cf. the old debate on the hard-sphere 
crystal 14 ) , presumably due to very small free-energy dif- 
ferences. For freely rotating ellipsoids a recent work has 
shown the existence of a different crystal packing consist- 
ing of a simple monoclinic lattice with a basis of two el- 
lipsoids with different orientation*^ (of course this struc- 
ture cannot occur in our parallel-particle model) . In our 
case of parallel PHSE particles with a < 1.5, simula- 
tions with between 8 to 12 layers along the z direction 
may give structures with any type of stacking depending 
on how the initial configuration is prepared (in the limit 
a = 1, we know that ABC stacking is preferred 14 ). How- 
ever, in the limit a — > 00, the AAA structure is clearly 
more stable, and it is expected that somewhere in the 
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interval 1.5 < a < oo there is a change in the nature of 
the crystal phase. Some preliminary simulations indicate 
that the critical value may be between 1.5 and 2.0 (e.g. a 
crystal with a = 1.5, when prepared with AAA stacking, 
shows a strong tendency to evolve towards random stack- 
ing, whereas the system with a — 2.0, when prepared 
with ABC stacking, evolves towards AAA stacking), but 
free-volume theory gives a rather higher value (see later) . 

Now the central question of the present paper is 
whether there exists an intermediate, stable fluid phase 
between the low-density nematic and the high-density 
crystal phase. In the limit a = 1 (parallel ellipsoids) 
ones knows for certain that there are no phases with par- 
tial (smectic or columnar) order—. By contrast, in the 
opposite limit a = oo (parallel cylinders) previous MC 
simulations^ indicate that the nematic phase changes to 
a smectic phase via a continuous transition. Clearly in 
our interaction model, which interpolates between these 
two limits, there must be a value of a beyond which 
the smectic phase becomes stable. Actually this is the 
case. In Fig. [2] cases where a smectic phase has been 
identified are indicated by a corresponding label (S for 
smectic). Evidence for the S phase comes from density 
distributions (see later), as the pressure shows no sign 
of a N-S transition, pointing to a continuous transition. 
We should note that recent simulations on freely-rotating 
hard spherocylinders in the limit of infinite aspect ratio 
have shown that the nematic-smectic transition is of first 
order—. Thus, we can conclude that the parallel align- 
ment constraint is responsible for the second-order nature 
of the N-S transition of hard cylinders, and this could 
also be the case in our model. By increasing the pressure 
further (compression run), the smectic phase transforms 
into a defected crystal phase via a first-order phase tran- 
sition. 

As mentioned above, our simulations indicate that a 
smectic phase stabilises for finite a, i.e. in the range 
a c < a < oo, where a c is some critical value. Indi- 
rect evidence comes from computation of density pro- 
files, smectic order parameter (not shown) and correla- 
tion functions. In Fig. |4] the evolution of the density 
profile p(z), as density is increased, is shown for the case 
a = 1.5. The onset of smectic order is clearly identified 
by the smooth appearance of density peaks. The density 
profiles exhibit a clear stratification at a packing fraction 
7y NS ~ 0.43. That the high-density phase in question is 
smectic, and not crystalline, can be concluded from the 
radial distribution functions along and perpendicular to 
the director, shown in Figs. [3](a) and (b), which point to 
in-plane fluid-like correlations in the intermediate phase. 

Analysis of the case a = 1.2, in contrast, suggests no 
evidence for a smectic phase: the nematic fluid freezes 
directly into a crystal (see Fig. [2]), similarly to the case 
of ellipsoids. The intermediate case a = 1.3, however, 
does show signs of smectic stability. >From all the in- 
formation collected for the various systems analysed, we 
estimate the critical value of a associated with the on- 
set of stability of the smectic phase to be in the range 
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Figure 4: Density profiles p(z) for the PHSE fluid with a 
1.5; (a) r\ = 0.39; (b) 77 = 0.43; and (c) 77 = 0.45. 



a c = 1.2-1.3. 

Fig. [5] summarises our results in the form of a phase 
diagram of packing fraction rj vs. inverse exponent a -1 . 
For a > 1.2 the smectic phase (S in the graph) is sta- 
ble. Simulation data for the limit a = 00 (cylinders, 
filled triangles and squares) are taken from simulations 
by Veerman and Frenkel—, and are essentially exact as 
they were inferred from simulations incorporating free- 
energy calculations. Data for hard spheres (filled circles) 
are taken from Hoover and Ree^. In the other cases 
transition densities for the nematic-smectic, nematic- 
crystal and smectic-crystal transitions are only approxi- 
mate. As already mentioned, the first (indicated by open 
circles) is of second order. The others (open triangles 
and squares) are of first-order, with a wide density gap; 
since no free energies were computed, we only plot, us- 
ing vertical bars, approximate limits of metastability of 
the two phases involved (the symbols are only the aver- 
age packing-fraction values inferred from the estimated 
limits of metastability). 

A second aim of our study is to rationalise the findings 
from computer simulation using theoretical models. As 
mentioned in Sec. Ill, we have considered three theoreti- 
cal models: a low-density virial expansion up to third or- 
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Figure 5: Phase diagram in the plane packing fraction 77- 
inverse shape parameter a' 1 . Continuous lines: coexistence 
boundaries for the smectic-solid transition from PL theory 
for the smectic and FV theory for the solid; shaded region: 
two-phase region from previous approximations; dashed line: 
nematic-smectic spinodal from PL theory; dotted: nematic- 
smectic spinodal from V3 theory; filled triangles: simula- 
tion results for coexistence packing fractions of smectic- 
solid transition in parallel cylinders (a = 00) from Veer- 
man and Frenkel^, and for liquid-solid coexistence in hard 
spheres^ (a = 1); filled square: nematic-smectic spinodal 
from computer simulation of Veerman and FrenkeiS'; open 
circles: our simulation estimates for nematic-smectic spin- 
odal; open squares: our simulation estimates for first-order 
smectic-solid transition; open triangles: our simulation esti- 
mates for first-order nematic-solid transition. In the latter 
two cases the vertical bars only indicate limits of metasta- 
bility. Labels indicate stable phases; N, nematic; S, smectic; 
ABC and AAA, crystalline solids with corresponding symme- 
tries. Vertical dot-dashed line: approximate limit of degen- 
eracy of ABC and AAA structures within FV theory. 



der in density for the low-density nematic phase (V3), a 
resummed virial expansion of the Parsons-Lee type (PL) 
for intermediate densities, and free- volume (FV) theory 
for the high-density crystal. The V3 theory was initially 
used to analyse possible bifurcations of the nematic phase 
to smectic or columnar phases^. The qualitative results 
as far as the continuous nematic-smectic transition is 
concerned are the same for both V3 and PL theories; 
they only differ quantitatively, as can be seen in Fig. [H 
For large values of a -1 (the only available from simula- 
tions), the MC data are bracketed by the two theories, 
but the trend that the packing fraction at bifurcation 
increases with a -1 in this regime is captured correctly. 
This increase is basically due to the decrease in particle 
volume as cylinders change to ellipsoids. The differences 
between the V3 and PL predictions can be traced back 
to their different treatment of correlations: V3 includes 
three-body correlations exactly, but the remaining terms 
are neglected, while PL only includes the exact two-body 



correlations but approximately resums the higher-order 
density correlations. 

The PL theory, together with the FV theory for the 
crystal, were used to compute the nematic-solid and 
smectic-solid transitions. These are first-order transi- 
tions, with a wide density gap. The two-phase region 
merges with the nematic-smectic spinodal coming from 
lower densities at a critical end-point located at a c ~ 1.2; 
this value approximately agrees with that inferred from 
the simulations. For lower values of a there is direct coex- 
istence between nematic and solid phases, and the theory 
also seems to agree with simulations in the limit of hard 
spheres (this agreement is fortuitous as it is well known 
that the PL theory is inappropriate to model spatial cor- 
relations in the hard-sphere fluid; the weighted-density 
theory-i^ or a theory based on FMT— are known to be 
more adequate). In the opposite limit (cylinders), how- 
ever, there is a big discrepancy with simulation, not only 
in the location of the nematic-smectic spinodal and the 
smectic-solid phase transition, but also in the density gap 
of the latter which is significantly overestimated as com- 
pared to the MC simulations. In this limit the weighted- 
density type theories, developed for a HSC fluid^, are 
obviously more appropriate to predict the spatial corre- 
lations. The PL and FV theories were also used to com- 
pute the smectic-columnar and columnar-solid transi- 
tions. These results are not plotted in the phase diagram 
since the region of columnar stability is always preempted 
by crystallization directly from the smectic phase; this 
feature also agrees with simulations which, as mentioned 
before, do not seem to give conclusive evidence for stable 
columnar ordering in this model. 
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Figure 6: Equation of state for the PHSE system with a = 
1.5. Symbols: MC results; continuous line: nematic branch as 
obtained from CS approximation; dashed line: smectic branch 
from PL theory; dotted line: solid branch according to FV 
theory. 

The quality of the different approximations, in the den- 
sity range where each of them was used, can be checked 
by examining the equation of state (EOS). This is done 
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in Fig. [6j which refers to the case a — 1.5. The fig- 
ure also contains the MC simulation data. The PL the- 
ory, which reduces to the Carnahan-Starling approxima- 
tion for the nematic phase (continuous line), represents 
correctly the EOS up to packing fractions of i] ~ 0.4, 
just before the transition to the smectic phase occurs 
(?7 NS ~ 0.43). The PL theory is not as accurate for the 
smectic phase (dashed line), since it overestimates the 
pressure. The nematic-smectic transition density (bifur- 
cation point) is reasonably well reproduced. In the solid 
branch (dotted line), the equation of state is accurately 
represented (in fact better as the density approaches the 
close-packing limit), as expected, but the smectic-solid 
transition is not correctly reproduced (value of pressure, 
location of transition densities and density gap), due to 
the defects in the smectic and solid equations of state. 



V. CONCLUSIONS 

We can conclude from our study that there is clear 
evidence that in the PHSE fluid smecticity begins to be 
favoured entropically beyond a critical value of the expo- 
nent a c ~ 1.2-1.3; this is surprisingly close to the value 
corresponding to ellipsoids (a = 1). Therefore, the fluid 
of parallel hard ellipsoids does not seem to be patholog- 
ical or special in not exhibiting a stable smectic phase; 
rather, there is a family of particle shapes beyond the 
ellipsoid that do not possess stable smectic phases, al- 
though the extent of the family in parameter space is 
small. 

The explanation for this behaviour lies, of course, in 
the packing efficiency of these hard-particle fluids. The 
formation of the smectic phase is the result of a deli- 
cate packing effect directly related to particle shape. In 
the HE fluid, the ellipsoidal geometry creates a large ten- 
dency for particles to interlock at their ends, which favors 
the stabilisation of the crystal phase with detriment to 
the smectic phase. This effect is not relevant in the case 
of HSC and sets in at some intermediate particle shape. 

Our simulation results seem to indicate that the colum- 
nar phase is not stable for this family of parallel rods. 
But free-energy computations will be needed to settle 
this question. However, in common with the case of par- 
allel cylinders, it is not likely that the columnar phase 
will be stable in our model. All these conclusions corre- 
spond to the model of parallel particles. Analysis should 
be extended to consider particles with free orientational 
degrees of freedom. Given the high degree of orienta- 
tional order of the smectic phases, we do no think the 
conclusions drawn from our study will be changed sub- 
stantially. 

As a byproduct of our analysis, we have assessed the 
validity of various theoretical approximations by compar- 
ison with the simulation results. A combination of free- 
volume and PL theories qualitatively predicts the correct 
phase stability, but the agreement is far from quantita- 
tive, except in the parameter region close to the spheres; 



in particular, the value of shape anisotropy beyond which 
smectic stability sets in is quite in agreement with sim- 
ulations. Use of the PL theory to consistently describe 
all phases is not adequate, since this theory progressively 
degrades as order builds up in the system: therefore, the 
agreement is probably fortuitous. 

An obvious avenue to improve the theoretical treat- 
ment is to use a more sophisticated theory, which is still 
to be formulated. A promising approach is FMT, which 
one hopes would correctly predict the relative stability of 
non-uniform phases; this belief is based on FMT calcula- 
tions applied to the particular case of parallel hard cylin- 
ders, as shown in Ref A However, the extension of this 
theory to superellipsoidal geometry would demand the 
use of strong approximations, in the line of Refi^» where 
only the first terms of the excess free-energy asymptotic 
expansion with respect to k _1 (with K the particle as- 
pect ratio) are taken into account. This procedure could 
disrupt the high predictive power of FMT as concerns 
the precise location of phase transitions and the relative 
stability of different non-uniform phases. Further inves- 
tigation of this problem would be needed but is left for 
future work. 

Interesting questions on the effect of particle shape on 
the depletion interaction between two anisotropic bod- 
ies mediated by small spherical particles can be studied 
within the present model. Recent studies have shown 
that particle geometry has a strong effect on the deple- 
tion forces^. The present model can be used to tune the 
geometry of the particle with a single parameter, from 
ellipsoids to cylinders, and allows the study of the evo- 
lution of the depletion potential with respect to a when 
two superellipsoids are immersed in a sea of small hard 
spheres. Work in this direction is in progress. 
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Appendix A: SPINODAL INSTABILITY OF N-S 
TRANSITION 



The spinodal instability coincides with the location of 
the continuous N-S phase transition. The instability con- 
dition can be calculated by solving the following set of 
equations 



1 - pc(r], q) = 0, 



dc(r), q) 



0. 



(Al) 
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where c(r), q) is the Fourier transform of the direct corre- 
lation function. The equations are to be solved for rf and 
q* , which are the value of the packing fraction and the 
wave number of the smectic phase, both at bifurcation. 
The second equation refers to the absolute minimum of 
—c(t], q) as this is a strongly oscillatory function of q. 

To calculate the spinodal curve a model for the direct 
correlation function is needed. The third virial approach, 
Eqn. (H]), gives 



c(»?,«)=/(g)+pA(g), 



(A2) 



where the Fourier transform of the Mayer function f(q) = 
J dre tqz f(r) gives the following result for a PHSE particle 
with unit breadth and height: 



/(«) 



4tt 
<1 



dzz sin 



9(1 



v 2a\!/(2a) 



and 



- A(g) = / dre^f(r)AV(r), (A4) 

AV(r) = J dr'/(r')/(r - r'). (A5) 

For each value of q (which is understood to be given in 
units of Co) , taken from a set of equally-spaced points, 
the function A(q) was calculated by MC integration for a 
fixed a. Then the equations IjAip were solved to find the 
spinodal curve rj*(a). Finally, the Parsons-Lee approach 
gives c(r),q) = ±* H s(?7)/(?)- 



Appendix B: PARAMETERIZATIONS OF THE 
PARSONS-LEE THEORY 



The density distributions for smectic and columnar 
symmetries were parameterized using two different ap- 
proaches. The first, based on a Fourier expansion, is 
more adequate for relatively low mean densities, while the 
other, a parameterization based on two parameters (one 
being a measure of the width of the density peaks, the 
other being the smectic or columnar lattice parameters - 
periods), is more appropriate for high-density phases, as 
the numerical convergence of the minimization schemes is 
much easier when the parameter space is reduced to two 
variables. We have used the following parameterization 
for the density profile of the smectic phase: 



p(z) 



Po 



exp (A cos qz), 



(Bl) 



where po = d" 1 J Q dzp(z) is the mean density, q = 2n/d 
the wave number (d being the smectic period), A the 
minimization parameter, and IqW the zeroth order mod- 
ified Bessel function. Inserting this expression into the 
Parsons-Lee functional, Eqn. ((Sj) , we find the following 



expression for ip cx = f^F^/N , the excess part of the free- 
energy per particle and unit thermal energy: 



2B^pl' 



dzp(z) 



dz'piz'mz-z 1 ), 

(B2) 



where B% = 4vo(a) (since the volume of the refer- 
ence HS particle is made to coincide with that of the 
PHSE). Also we have defined f(z) = J dr±f(r±) as the 
integrated Mayer function over the transverse area, with 
r_L = (x, y). Note that the HS free-energy per particle is 
^nsiv) = (4 — 3/7)77/ (1 — rf) 2 . By calculating the integral 



involved in Eq. (|B2[) we find explicitly 



(A3) where 



2 \a'2a )' 



(B3) 



(B4) 



X (A;a) = I - 2 (A) J dz(l - z 2a )^ a I (^2Acos l± 

(B5) 

and q* = 2-na^jd. Finally, the ideal part of the free- 
energy per particle and unit thermal energy, tp\& = 
PTid/N, can be found as 



ipid = In 7/ — 1 + A 



h (A) 
*o(A) 



In J (A), 



(B6) 



where I\{x) is the first order modified Bessel function. 
We have minimized the total energy ip — ipid + Vox with 
respect to A and d for a fixed value of r\. 

For the columnar phase, the parameterization chosen 
is a sum of Gaussian peaks centered at the sites of the 
triangular lattice: 



AoPo^ccll 



^exp 

k 



-A (rj_ - R k )" 



(B7) 



and normalized in such a way that integration over the 
unit cell of area ^4 C oii = V3a 2 /2 (a being the lattice pa- 
rameter of the triangular cell) gives the mean density po- 
The position of the lattice sites are Rk = fciai + &2a2 

(ki G Z) with a„ = ^ (— 1)™J being the vectors of 

the triangular lattice. The Gaussian width is controlled 
by the parameter Ao which, together with a, define the 
set of minimization variables. The excess part of the 
free-energy per particle can be calculated as 



</?cx 





1 - (ra*) 2Q ] 


x / drr 


Jo 





,l/(2o) 



(B8) 
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where 7 = Xqo 2 , a* — a/ao, i?£ = |Rk|/a, and 

^ccii = ( a * ) (the un it cell area in dimensionless 
units) . Correspondingly the ideal part of the free-energy 
per particle is defined by 



In 



V3rn/(2ir) 



1 



27T pOO 

d<f) I drr 
n Jo 



x e r In 



E 1 



(B9) 



An useful approximation for the ideal part when 7 > 1 
can be obtained by taking RJj. — > in Eqn. (|B9j) . with 
the result 



ip id = In 



V3r n /(2n) 



(BIO) 



Also, the excess part (|B8|) can be approximated by taking 
only the terms k = (0, ±1), k = (±1, 0), k = (1,-1) and 
k = (—1,1) in the sum (i.e. considering only the nearest 
neighbours of a given site) . 

Now we proceed to describe the other parametriza- 
tion used for minimization at low densities: a truncated 
Fourier expansion of the density the density profile. This 
expansion for the smectic symmetry is 



p(z) = poip(z) = po 



N 



1 + 2_. Pn cos(gnz) 



n>0 



(Bll) 



with p n the Fourier amplitudes. Using this parameter- 
ization, the excess part of the free-energy per particle 
reads 



A* 



Ye 



(B12) 



«>o 



where 9 r , 



(1 + S n<0 )/2 and 



T n = 



4C Q 



dzz sin 



(1- 



2a\V(2a) 



(B13) 



In Jo 

with q n — q<7 n, while the ideal part is 

tp id =lnry- 1 + 2 / dzil>(dz)1nip(dz). (B14) 
Jo 



For the columnar symmetry the Fourier expansion of 
the density distribution reads 



p(r x ) = Poip(r x ) = pa 



Ni,N 2 

1 + Pnxn 2 
ni : U2>0 



'2miix\ ( lixn-iy 
x cos [ =r- cos ' 



a\/l\ 



(B15) 



where p ni n 2 are the two-dimensional Fourier amplitudes. 
After insertion of Eqn. (|B15|) into the Parsons-Lee ex- 
pression for the excess free energy, we obtain 

fcx^^miv) E °r lin2 T ni n 2 pl in2 , (B16) 

nl,n2>0 

for the excess free-energy per particle, where 6 ni „ 2 = 

(^ni+na.O + 8n u + ^n 2 ,0 + l)/4, and 
»1 

T nin2 = 4C Q / dzz (1 - z 2 «) 1/(2Q) J (g nin2 z),(B17) 
Jo 

with Jq(x) the zeroth order Bessel function and 

9nin 2 = — \ 1" n i- The ideal part of the free-energy 

a* V 3 

per particle for columnar symmetry can be calculated as 
tpid = In f] - 1 



dx dy 
OVott Jo Jo 



x ip(x,y)lntj}(x,y). 



(B18) 



This time the total free-energy per particle should be 
minimized with respect to the lattice parameter a and 
also with respect to the Fourier amplitudes p ni n 2 - While 
the Fourier parameterizations were used for the min- 
imization of the total free-energy of the smectic and 
columnar phases at low densities (up to 77 ~ 0.55), the 
corresponding parameterizations using only two param- 
eters, described above, were used to minimize the free 
energies per particle of these phases at higher densities. 
Both approaches were checked for consistency at inter- 
mediate densities. 
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